The diprate package is available from GitHub here: dipDRC
library(diprate)
options(stringsAsFactors=FALSE)
static <- read.csv("../data/from_CW_via_Slack_20210507/StaticDF.csv", row.names=1)
static <- static[order(static$Culture_Type,static$Cell_Line),]
static[static$Cell_Conc==0,"Cell_Conc"] <- 1
cell_lines <- unique(static$Cell_Line)
static <- static[!(static$Cell_Line=="CORL279" & static$RLU < 10),]
par(mfrow=c(2,5))
linear_models <- lapply(cell_lines, function(cl) {
dat <- static[static$Cell_Line==cl,]
culture_type <- unique(dat$Culture_Type)
m <- lm(log2(RLU) ~ log2(Cell_Conc), data = dat)
plot(log2(RLU) ~ log2(Cell_Conc), data=dat, main=paste0(cl," (",culture_type,")"), xlim=c(0,14), ylim=c(7,18))
abline(m, col="blue")
return(m)
})
names(linear_models) <- cell_lines
Assuming first part of data is at some minimum value and not associated with any actual cell count (lower limit of detection), which would result in slope=0 (values are constant until some minimum number of cells is achieved).
lagLine <- function (x, lower, slope, br = 64) sapply(x, function(z) ifelse(z <= br, lower, lower + (z-br) * slope))
fitLagLin <- function(x, y, start_list=list(lower=5, slope=1, br=1))
nls(y ~ lagLine(x=x, lower, slope, br),
start = start_list,
algorithm="port",
control=nls.control(maxiter=500)
)
par(mfrow=c(2,5))
lag_linear_models <- lapply(cell_lines, function(cl) {
m <- tryCatch({
dat <- static[static$Cell_Line==cl,]
culture_type <- unique(dat$Culture_Type)
m <- fitLagLin(log2(dat$Cell_Conc), log2(dat$RLU))
}, error=function(e) { return(e) })
# if(culture_type == "Adherent")
# {
# xr <- c(4,12)
# } else {
# xr <- c(6,14)
# }
xr <- c(0,14)
plot(log2(RLU) ~ log2(Cell_Conc), data=dat, main=paste0(cl," (",culture_type,")"), xlim=xr, ylim=c(7,18))
if(class(m)[1] != "nls")
{
m <- lm(log2(RLU) ~ log2(Cell_Conc), data=dat[dat$Cell_Conc>1,])
abline(m, col="blue", lwd=2)
} else {
curve(from=0.5,to=18, lagLine(x, lower=coef(m)['lower'],
slope=coef(m)['slope'],
br=coef(m)['br']),
col="blue", lwd=2, add=TRUE)
}
# text(12, 9, paste("Adj R2 ="))
return(m)
})
names(lag_linear_models) <- cell_lines
Assume all luminescence values with cells produce detectable signal.
par(mfrow=c(2,5))
linear_models <- lapply(cell_lines, function(cl) {
dat <- static[static$Cell_Line==cl & static$Cell_Conc >1,]
culture_type <- unique(dat$Culture_Type)
m <- lm(log2(RLU) ~ log2(Cell_Conc), data = dat)
plot(log2(RLU) ~ log2(Cell_Conc), data=dat, main=paste0(cl," (",culture_type,")"), xlim=c(0,14), ylim=c(7,18))
abline(m, col="blue")
return(m)
})
names(linear_models) <- cell_lines
Assuming the lowest number of cells is above the threshold of detection, will remove the no cells control and fit remaining data.
dat <- static[static$Cell_Conc > 1,]
m2 <- lme4::lmList(log2(RLU) ~ log2(Cell_Conc) | Cell_Line, data=dat)
Registered S3 methods overwritten by 'lme4':
method from
cooks.distance.influence.merMod car
influence.merMod car
dfbeta.influence.merMod car
dfbetas.influence.merMod car
f2 <- coef(m2)
r2 <- unlist(summary(m2)$adj.r.squared)
par(mfrow=c(2,5))
temp <- lapply(cell_lines, function(cl) {
dtp <- dat[dat$Cell_Line==cl,]
culture_type <- unique(dtp[dtp$Cell_Line==cl,'Culture_Type'])
if(culture_type == "Adherent")
{
xr <- c(5,12)
} else {
xr <- c(7,14)
}
plot(log2(RLU) ~ log2(Cell_Conc),
data=dtp,
main=paste0(cl," (",culture_type,")"),
xlim=xr, ylim=c(7,18))
abline(m2[[cl]], col="blue", lwd=2)
text(xr[1]+0.25, 17, pos=4, paste("slope =",signif(f2[cl,2],3)))
text(xr[1]+0.25, 16, pos=4, expression(R^2))
text(xr[1]+1, 16, pos=4, paste("=", signif(r2[cl],3)))
})
fitLagLin1 <- function(x, y, start_list=list(lower=5, br=1))
nls(y ~ lagLine(x=x, lower, slope=1, br),
start = start_list,
algorithm="port",
control=nls.control(maxiter=500)
)
par(mfrow=c(2,5))
lag_linear1_models <- lapply(cell_lines, function(cl) {
m <- tryCatch({
dat <- static[static$Cell_Line==cl,]
culture_type <- unique(dat$Culture_Type)
m <- fitLagLin1(log2(dat$Cell_Conc), log2(dat$RLU))
}, error=function(e) { return(e) })
# if(culture_type == "Adherent")
# {
# xr <- c(4,12)
# } else {
# xr <- c(6,14)
# }
xr <- c(0,14)
plot(log2(RLU) ~ log2(Cell_Conc), data=dat, main=paste0(cl," (",culture_type,")"), xlim=xr, ylim=c(7,18))
if(class(m)[1] != "nls")
{
m <- lm(log2(RLU) ~ log2(Cell_Conc), data=dat[dat$Cell_Conc>1,])
abline(m, col="blue", lwd=2)
} else {
curve(from=0.5,to=18, lagLine(x, lower=coef(m)['lower'],
slope=1,
br=coef(m)['br']),
col="blue", lwd=2, add=TRUE)
}
# text(12, 9, paste("Adj R2 ="))
return(m)
})
names(lag_linear1_models) <- cell_lines
Cells seeded on Dec 15 Drug added:2 pm Dec 16; 60 + 20 µl of drug Data acquisition began:
lcc <- read.csv('../data/from_CW_via_Slack_20210507/20201216_Lum_CellCounts_Thunor.csv', row.names=1, as.is=TRUE)
lcc <- lcc[,-1]
lcc$uid <- paste(lcc$upid,lcc$well,sep="_")
lcc <- lcc[order(lcc$uid,lcc$time),]
lcc <- lcc[lcc$time <= 96,]
lcc_cell_lines <- unique(lcc$cell.line)
ctrls <- lapply(lcc_cell_lines, function(cl) lcc[lcc$cell.line==cl & lcc$drug1.conc==0,])
names(ctrls) <- lcc_cell_lines
par(mfrow=c(2,2))
invisible(lapply(names(ctrls), function(n) do.call(plotGC,
append(getGCargs(ctrls[[n]], dat.col=c("time","Cell_Count","uid")),list(main=n, leg=FALSE)))))
par(mfrow=c(2,2))
invisible(lapply(names(ctrls), function(n) do.call(plotGC,
append(getGCargs(ctrls[[n]], dat.col=c("time","RLU","uid")),list(main=n, leg=FALSE)))))
sumc <- do.call(rbind, lapply(lcc_cell_lines, function(cl)
{
counts <- sapply(unique(ctrls[[cl]]$time), function(i)
sum(ctrls[[cl]][ctrls[[cl]]$time==i,"Cell_Count"]))
times <- unique(ctrls[[cl]]$time)
data.frame(cell.line=cl, time=times, cell.count=counts)
}))
sumc$popdubs <- log2norm(sumc$cell.count,sumc$cell.line)
sumc$cell.line <- as.character(sumc$cell.line)
# h1048_sumc <- data.frame(time=unique(h1048$time), cell.count=h1048_sumc)
# plot(log2(cell.count)-log2(cell.count)[1] ~ time, data=h1048_sumc, type="l", ylab="Population doublings")
par(mfrow=c(3,2))
invisible(lapply(lcc_cell_lines[lcc_cell_lines!="WM88"], function(cl)
{
dtp <- ctrls[[cl]]
invisible(do.call(plotGC, append(getGCargs(dtp, dat.col=c("time","Cell_Count","uid")),
list(main=paste0(cl,", cell count"), leg=FALSE))))
invisible(do.call(plotGC, append(getGCargs(dtp, dat.col=c("time","RLU","uid")),
list(main=paste0(cl,", lum"), leg=FALSE, ylim=c(-1,3)))))
lines(sumc[sumc$cell.line==cl,"time"], sumc[sumc$cell.line==cl,"popdubs"], lwd=3)
}))
lumo <- read.csv("../data/from_CW_via_Slack_20210507/050819_RTglowDF.csv")
lumo <- lumo[,-1]
lumo_cell_lines <- unique(lumo$cell.line)
lumo_ctrl_dat <- lapply(lumo_cell_lines, function(cl) lumo[lumo$cell.line==cl & lumo$drug1.conc==0 & lumo$drug1!="DMSO",])
names(lumo_ctrl_dat) <- lumo_cell_lines
par(mfrow=c(2,4))
invisible(lapply(lumo_cell_lines, function(cl)
{
dtp <- lumo_ctrl_dat[[cl]]
# plot(dtp$TotHour, dtp$RLU)
invisible(do.call(plotGC, append(getGCargs(dtp, dat.col=c("TotHour","RLU","well"), arg.name = c("time", "cell.count", "ids")),
list(main=cl, leg=FALSE))))
}))
corl279 <- lumo_ctrl_dat[['CORL279']]
a <- subset(corl279, grepl("12", corl279$well) )
invisible(do.call(plotGC, append(getGCargs(a, dat.col=c("TotHour","RLU","well"), arg.name = c("time", "cell.count", "ids")),
list(main="CORL279", leg=FALSE))))
dms53 <- ctrls[["DMS53"]]
sapply(lag_linear1_models, function(x) coef(x)[1]-coef(x)[2])
DMS114.lower DMS454.lower DMS53.lower H1048.lower H841.lower CORL279.lower H524.lower H526.lower
5.0471261 6.0059948 7.4546528 6.6395654 4.7169988 2.6679066 -0.4606302 1.1555825
par(mfrow=c(1,3))
invisible(lapply(lcc_cell_lines[lcc_cell_lines!="WM88"], function(cl)
{
dtp <- ctrls[[cl]]
invisible(do.call(plotGC, append(getGCargs(dtp, dat.col=c("time","RLU","uid")),
list(main=paste0(cl,", lum"), leg=FALSE, ylim=c(-1,3)))))
lines(sumc[sumc$cell.line==cl,"time"], sumc[sumc$cell.line==cl,"popdubs"], lwd=3)
}))
d <- read.csv("../data/20201216_Lum_Image_Run/20201216_Lum_CellCounts_TOTAL.csv", row.names = 1)
d$uid <- paste(d$Plate_Name,d$Well,sep="_")
d <- d[order(d$uid,d$TotHour_Lum),]
d <- d[d$TotHour_Lum <= 96,]
par(mfrow=c(1,3))
invisible(lapply(lcc_cell_lines[lcc_cell_lines!="WM88"], function(cl)
{
dtp <- d[d$Cell_Line==cl & d$Drug=="AZD-1152", ]
invisible(do.call(plotGC, append(getGCargs(dtp, arg.name = c("time", "cell.count", "ids", "rep"), dat.col=c("TotHour_Lum","RLU","uid","Drug_Conc")),
list(main=paste0(cl,", lum"), ylim=c(-1,3)))))
}))